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Auxiliary-field  quantum  Monte  Carlo  method  for  strongly  paired  fermions 
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We  solve  the  zero-temperature  unitary  Fermi  gas  problem  by  incorporating  a  BCS  importance  function  into 
the  auxiliary-held  quantum  Monte  Carlo  method.  We  demonstrate  that  this  method  does  not  suffer  from  a  sign 
problem  and  that  it  increases  the  efficiency  of  standard  techniques  by  many  orders  of  magnitude  for  strongly 
paired  fermions.  We  calculate  the  ground-state  energies  exactly  for  unpolarized  systems  with  up  to  66  particles 
on  lattices  of  up  to  273  sites,  obtaining  an  accurate  result  for  the  universal  parameter  § .  We  also  obtain  results  for 
interactions  with  different  effective  ranges  and  find  that  the  energy  is  consistent  with  a  universal  linear  dependence 
on  the  product  of  the  Fermi  momentum  and  the  effective  range.  This  method  will  have  many  applications  in 
superfluid  cold  atom  systems  and  in  both  electronic  and  nuclear  structures  where  pairing  is  important. 

DOI:  10.1 103/PhysRevA.84.061602  PACS  number(s):  67.85.Lm,  03.75.Hh,  03.75. Ss,  24.10.Cn 


The  study  of  strongly  interacting  Fermi  systems  is  one 
of  the  central  themes  and  major  challenges  in  physics. 
Superfluidity  in  unpolarized  cold  atomic  Fermi  gases,  which 
has  been  demonstrated  both  experimentally  and  theoretically, 
provides  a  prototypical  example.  The  experimental  ability 
to  use  a  Feshbach  resonance  to  adjust  the  strength  of  the 
potential  between  the  atoms  allows  an  exploration  of  the 
physics  over  many  length  scales.  A  particularly  interesting 
regime  is  at  unitarity  where  the  scattering  length  diverges  and 
the  effective  range  of  the  potential  is  very  small  compared  to 
the  interparticle  spacing.  Since  the  particle  density  provides  the 
only  length  scale,  the  ground-state  energy  Eq  is  proportional 
to  the  free  Fermi  gas  energy  EFq, 

Eo  =  %Efg  ■  (1) 

The  ability  to  quantitatively  understand  the  properties  of 
this  system  symbolizes  a  great  triumph  of  many-body  physics. 
Many  experiments  and  calculations  have  been  performed  for 
the  unitary  Fermi  gas.  Initial  qualitative  agreement  was  found 
between  theoretical  [1,2]  and  experimental  [3-5]  results.  More 
precise  recent  experiments  have  yielded  q  =  0.39(2)  [6]  and 
0.41(1)  [7],  with  smaller  values  obtained  very  recently  by 
Zwierlein  et  al.  [8],  Fixed-node  diffusion  Monte  Carlo  (DMC) 
calculations  [1,2,9-12]  have  always  included  a  Bardeen- 
Cooper-Schrieffer  [13]  (BCS)  trial  wave  function  to  guide  the 
Monte  Carlo  walk  and  provide  the  fixed  node  constraint  [14] 
needed  to  overcome  the  fermion  sign  problem.  As  is  well 
known,  these  calculations  provide  an  upper  bound,  with  the 
current  best  value  f  =  0.383(1)  [11,12]. 

In  this  paper,  we  show  that  exact  calculations  can  be 
performed  to  accurately  determine  the  ground-state  properties 
of  the  unpolarized  Fermi  gas.  Furthermore,  we  perform 
calculations  with  small  but  finite  effective  range  and  find 
results  consistent  with  a  universal  linear  dependence  of  the 
energy  upon  the  Fermi  momentum  times  the  effective  range: 
E/EFq  —  §  +  Skpre  +  •  •  • . 

A  method  is  introduced  to  allow  the  use  of  a  BCS  trial 
wave  function  in  the  auxiliary-field  quantum  Monte  Carlo 
(AFQMC)  approaches  of  Zhang  and  coworkers  [15,16].  Using 
this  approach,  we  perform  calculations  with  several  forms  of 


the  kinetic  energy  term  that  all  give  the  correct  continuum 
limit  but  with  different  finite  effective  ranges  to  study  the 
convergence  with  particle  number  and  lattice  sizes  and  to 
obtain  the  dependence  of  if  on  the  effective  range.  An  exact 
result  is  obtained  for  the  value  of  £,  as  well  as  results  for  the 
slope  S. 

Quantum  Monte  Carlo  simulations  play  a  key  role  in 
addressing  the  challenge  of  strongly  interacting  Fermi  systems. 
The  AFQMC  method  has  been  applied  to  a  variety  of  systems 
in  several  fields.  With  equal  numbers  and  masses  of  up-  and 
down-spin  fermions  and  an  attractive  interaction,  there  is  no 
fermion  sign  problem.  The  formalism  developed  here  allows 
the  use  of  a  BCS  importance  function,  which  drastically  im¬ 
proves  the  efficiency  in  this  situation.  In  general  applications, 
a  sign  or  phase  problem  is  present,  which  is  controlled  by  a 
constraint,  also  using  the  importance  function  [15,16].  Hartree- 
Fock  or  free  Fermi  gas  (FG)  types  of  importance  functions  have 
been  used.  This  approach  has  been  shown  to  be  very  accurate 
in  many  condensed- matter  models  and  optical  lattices  [17], 
quantum  chemistry  [18],  and  solid-state  materials  [19].  The 
present  method  will  allow  the  use  of  BCS  importance  functions 
[or  antisymmetrized  geminal  power  (AGP)  in  chemistry]  to 
significantly  improve  our  ability  to  deal  with  the  sign  problem 
in  systems  where  pairing  is  important  and  enhance  the  capabil¬ 
ities  for  quantum  simulations  in  strongly  correlated  systems  in 
general. 

The  AFQMC  method,  in  both  zero-  [20]  and  finite- 
temperature  formulations  [21],  has  also  been  applied  to  the 
unitary  Fermi  gas.  Precise  results  require  simultaneously 
large  lattices,  so  the  system  is  dilute,  and  a  large  number 
of  particles  for  an  accurate  approach  to  the  continuum  and 
thermodynamic  limits.  Many  such  calculations  have  been 
performed  [20,22-25],  but  the  variance  of  the  method  limited 
the  results  to  a  relatively  small  number  of  particles  and  lattice 
sizes  so  that,  as  we  demonstrate  below,  the  results  are  unlikely 
to  have  converged  to  the  thermodynamic  limit. 

The  range  of  the  van  der  Waals  interaction  in  cold  atoms  is 
small  (e.g.,  about  3  nm  in  6Li  [26])  compared  to  interparticle 
spacing  so  that  the  short-range  structure  of  the  interaction 
is  unimportant;  the  results  are  completely  determined  by  the 
form  of  the  kinetic  energy  and  the  scattering  length.  For  an  N l 
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lattice,  the  equivalent  Hamiltonian  is 
H  =  ^3  H  if)s^mseke,k<ri~rm)  +  U  ^\i;t774  .  (2) 

^  k.j.m.s  i 


Here  \[rjs  is  the  destruction  operator  for  a  fermion  of  spin 
s  on  lattice  site  at  position  r j.  For  odd  lattice  sizes  Nk  = 
2NC  +  1 ,  the  k  are  given  by  =j-(nxx  +  nyy  +  nzz)  with 
— Nc  ^  nx,ny,nz  ^  Nc.  The  k  space  destruction  operators 
are  cks  =  e~,k'Tt\ jfjs,  and  the  density  operators  are 

nis  =  Wisfis- 

To  reach  the  continuum  limit,  we  need  to  take  the  limit  of 
zero  particle  density,  p  =  N/Nk  0,  in  the  context  of  the 
Hubbard  model  (i.e.,  replace  L  by  Nk).  Equivalently,  because 
of  scale  invariance,  we  can  think  of  the  system  as  a  discretized 
representation  of  a  supercell  with  fixed  size  L  and  take  the  k- 
space  cutoff  to  infinity.  In  either  case,  we  then  take  the  number 
of  particles,  N,  to  infinity.  In  this  limit,  only  the  behavior  of 
ek  for  k  <5jC  ^  is  important,  where  a  =  L/Nk  is  the  lattice 
spacing.  Thus,  a  variety  of  kinetic  energy  forms  can  be  used  as 
long  as  they  are  quadratic  for  k  much  smaller  than  the  cutoff. 
In  this  work,  we  present  results  for 


„(2) 


2  lA 


hlk 


ck 


2  7,2 


Trk 


2m 


pW], 


K  2m 

(h)  /i- 

e.  =  - ,  [3  —  cos (kxa)  —  cos(kva)  —  cos(k,a)] . 


(3) 


The  superscripts  2  and  4  indicate  the  highest  power  of  k,  while 
€<h>  is  the  Hubbard  model  hopping  kinetic  energy  offset  by  a 
constant  so  that  it  is  zero  at  k  —  0. 

For  two  particles,  the  Hamiltonian  is  separable,  and  the 
solution  of  the  Lippmann-Schwinger  equation  for  low-energy 
s-wave  scattering  gives  the  phase-shift  equation, 


where  P  indicates  the  principal  parts  integration  and  the  k 
space  sums  are  cut  off  by  the  lattice  spacing  a.  The  effective 
range  expansion  is 


kcot<5o  =  —a  1  +  \k2re  +  •  •  •  ,  (5) 

where  a  is  the  scattering  length  and  re  is  the  effective  range. 
Since  we  are  interested  in  the  unitary  limit,  we  adjust  U  to  have 
a  ~ 1  =  0.  Both  e[2!  and  e^n  have  nonzero  effective  ranges.  The 
extra  parameter  ft  in  c[4>  can  be  adjusted  to  make  the  effective 
range  zero.  The  values  for  the  parameters  are  given  in  Table  I. 

The  AFQMC  algorithm  uses  branching  random  walks  to 
project  the  ground  state  from  an  initial  trial  state  with  the 
imaginary-time  operator  exp[— Ht],  Because  the  interaction  is 
attractive,  there  is  no  fermion  sign  problem  for  equal  numbers 
of  up  and  down  fermions  studied  here,  and  no  path  constraint 
is  required.  A  walker  is  a  set  of  N  single-particle  orbitals. 
Initially,  the  orbitals  for  the  up-spin  particles  are  taken  to  be 
identical  to  those  for  the  down-spin  particles.  The  two-body 
interaction  term  is  broken  up  using  a  Hubbard-Stratonovich 
(HS)  transformation,  which  has  only  positive  weights  and 


TABLE  I.  The  parameters  that  give  infinite  scattering  length  for 
two  particles  in  an  infinite  lattice  for  the  various  kinetic  energies. 
The  fi  value  for  the  kinetic  energy  has  been  adjusted  to  give  zero 
effective  range,  re. 


Energy 

j  j  2 ma^ 

U  Ti2 

P 

rea  1 

4h> 

-7.91355 

-0.30572 

42' 

-10.28871 

0.33687 

44) 

-8.66605 

0.16137 

0.00000 

treats  the  up-  and  down-spin  particles  identically.  Therefore, 
the  up-spin  orbitals  remain  identical  to  the  down-spin  orbitals 
during  the  propagation.  We  show  below  that  the  usual  form  for 
a  singlet  paired  BCS  trial  function  also  gives  no  fermion  sign 
problem. 

The  walker  states  are  given  by  specifying  the  orbital 
coefficients.  These  can  be  specified  on  the  real  space  lattice 
<j)nj  or  as  momentum  space  coefficients  *  related  to  each 
other  by  a  discrete  Fourier  transform.  If  we  begin  with  real 
orbitals  on  the  real  space  lattice,  the  orbitals  remain  real  when 
propagated.  The  momentum  space  orbitals  therefore  satisfy 
4>n-k  =  4>„  k-  The  orbitals  are  orthonormalized  periodically. 
This  is  needed  to  limit  roundoff  error,  but  the  mathematical 
expressions  are  correct  without  it.  The  orthonormal  orbitals 
therefore  satisfy  &n  k&m- k  =  anc*  t^le  corresponding 

operators,  wns  —  /L/:  satisfy  {w —  8nm8SS'. 

The  walker  state  is 


N/2 

n= 1 


Because  N  N%  and  because  the  imaginary-time  history  of 
the  walk  need  not  be  retained,  this  formalism  is  much  more 
efficient  than  the  usual  lattice  formulations  for  the  ground  state 
of  dilute  gases. 

Using  a  discrete  HS  transformation  [27],  the  potential 
energy  propagator  is 


e~uJ2ini-tnnAt 


1 

2^ 


Gy({er},  Ar), 

M=±l 


Gy  ({a},  At)  =  exp 


(?) 

UOi  , 


where  tanh2  u  —  —  tanh(TAl).  The  kinetic  energy  propagator 
is 


GrlAr)  =  exp 


J]e*(4tc*t  +  4|cn)Af 

k 


(8) 


Given  a  choice  of  one  of  the  Nk  set  of  fields,  the  application 
of  the  Trotter  breakup  of  one  term  of  the  propagator  on  a  walker 
|  W)  gives  another  walker  |  W')  times  a  weight  w'({(j},W)  that 
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depends  on  the  set  of  HS  variables  {er}  and  |  W), 

G{{a},At)\W)  =  Gt  (y)  Gv({a},At)GT  (^W> 

w'({a},W)\W').  (9) 

The  propagation  above  consists  of  these  steps:  (1)  Multiply 
each  (pn  ji  by  exp(—  At).  (2)  Fast  Fourier  transform  them  to 
obtain  in  real  space.  (3)  Multiply  each  <j)n  j  by  exp(2«a,  — 
jU At).  (4)  Fast  Fourier  transform  them  to  obtain  <j)n  ^  in 
momentum  space.  (5)  Multiply  each  <j>n.k  by  exp(— ^ejAt). 
(6)  Update  the  weight  from  nonoperator  terms. 

Including  importance  sampling  reduces  the  fluctuations 
by  changing  the  sampling  so  that  it  is  nonuniform,  without 
biasing  the  results.  We  want  to  sample  walkers  |  W)  from 
(4't’I  W)(W\ijs(t)),  where 

+  At))  =  e~(H~ET)A,W(t)).  (10) 


The  contribution  of  a  walker  |  W )  at  the  next  time  step  is  then 


E 

{cr}=±l 


1  (ttr|G({q},AQ|W)  £rA,' 
2Nl  OPrlW) 


x 


1 

w({cr},  W) 


G({a], Af)|W). 


(ID 


We  want  to  sample  the  set  of  HS  variables  {a}  from  the 
unnormalized  probability  distribution  given  by  the  square 
brackets.  The  normalization  to  order  At2  is  the  local  energy 
expression  e~('^EL(W)+EL^w'^~Er'>A'  and  will  give  the  weight 
of  the  sampled  walkers.  Once  we  have  sampled  the  {cr }  values, 
the  new  walker  is  given  by  the  second  term  of  Eq.  (11).  We 
make  sure  that  regions  where  the  trial  function  are  small  are 
sampled  adequately  to  eliminate  trial-function  bias. 

The  particle  projected  BCS  state  is 


|BCS)  = 


(12) 


where  fa  =  /u /(  in  the  usual  notation.  The  overlap  of  a 
walker  with  the  BCS  state  is 

N/2  /  \  N/2 

(W|BCS>  =  (0|  f]  WniWn t  1°)  <13) 

n= 1  \  k  / 


The  contraction  needed  is  uinj/cl  =  4>„^ss'  ■  The  two  creation 
operators  in  the  BCS  pair  must  be  contracted  with  some  wm- j- 
and  some  w, 4,  giving  a  term 

Km  =  Wnjvmf  E ^4,^  (14) 

k 

^  '  &  n.—kfk&mk  ^  ’ 

k  k 


Taking  all  the  different  possible  full  contractions. 


(W|BCS)  =  detA  ,  (15) 

where  the  elements  of  the  4  x  t  matrix  A  are  the  Amn  of 
Eq.  (14). 

The  overlap,  Eq.  (15),  is  positive  when,  as  in  the  standard 
singlet  paired  BCS  solutions,  fa  ^  0.  We  write  detA  = 


det[BZ?t],  where  fit  js  the  Hermitian  conjugate  matrix  of 
B  and  the  matrix  elements  of  the  4  x  N%  matrix  B  are 
Bnk  —  (pn.k-y/Jk-  By  applying  the  Cauchy-Binet  theorem,  each 
of  the  determinants  of  the  submatrices  of  B  is  multiplied 
by  the  determinant  of  the  corresponding  Hermitian  conjugate 
submatrix.  The  determinant  of  A  is  a  sum  of  positive  terms,  so 
that  our  BCS  trial  function  gives  no  sign  problem. 

Two  estimates  of  the  energy  are  used:  The  growth  energy 
just  measures  the  growth  of  the  weights  in  the  random  walk. 
The  local  energy  can  be  calculated  using  contractions  similar  to 
those  above.  Other  observables  can  be  calculated  similarly.  A 
complete  derivation  for  a  general  BCS  state  will  be  published 
elsewhere.  Here  we  give  the  result 


(WjfflBCS)  , 

El(W)  =  - -  =  2tr[A~1C] 

(WjBCS) 

+  U  ^MA-^tejMA-'Et^)] 


-  tr ■[A~1E(q)A-lE\q)]  +  tr  [A~lD(q)]},  (16) 


where  Dnm)q)  —  ^ h .k  vq  fk^m.k-yq ■>  Cnm  — 

fk&m,k>  and  Enl)l  ({/  )  =  'Kk  tf\t.k  Yq  .k  ■  The 

matrix  elements  of  D  and  E  are  convolutions  which  are 
efficiently  computed  with  fast  Fourier  transforms.  The 
computational  cost  of  using  the  BCS  |4/j-)  is  similar  to  that 
using  a  single  Slater  determinant. 

We  have  calculated  the  ground- state  energy  for  different 
particle  numbers  and  lattice  sizes.  The  time-step  errors  have 
been  extrapolated  to  zero  within  statistical  errors,  and  walker 
population  biases  have  been  checked  and  were  found  to  be 
negligible  for  the  population  sizes  used.  The  imaginary  time 
step  is  «0. 01-0. 05  EJ 1 ,  the  total  propagation  time  is  of  order 
10-30  Epl,  and  2  000-20  000  random  walkers  are  used  in  the 
simulations. 

For  N  —  4,  we  found  that  the  use  of  BCS  importance 
functions  reduced  the  statistical  error  by  a  factor  of  10,  or 
lOOx  reduction  in  computer  time,  compared  to  the  usual  FG 
importance  function.  The  improvement  increased  to  1500x  for 
N  =  38  in  a  123  lattice.  For  larger  systems,  the  discrepancy 
is  much  larger  still;  indeed  the  statistical  fluctuations  from  the 
latter  are  such  that  often  meaningful  results  cannot  be  obtained 
with  the  run  configurations  described  above. 

In  Fig.  1,  we  summarize  our  calculations  of  the  energy  as 
a  function  of  p1^  where  p  =  N/Nj?,  and  the  particle  number 
is  N  —  38,  48,  or  66.  We  plot  §,  Eq.  (1),  where  we  have  in 

ft2 It2 

all  cases  used  the  infinite  system  free-gas  energy  EFq  —  | 
with  k\  —  3tt2-^-t  as  the  reference. 

DMC  calculations  have  found  converged  results  when 
using  66  particles  [11,12],  and  our  results  confirm  this. 
The  differences  between  38  and  66  particles  are  rather 
small.  Our  calculations  with  14  particles  show  a  significant 
size  dependence,  and  with  26  particles  the  effects  are  still 
noticeable.  These  are  not  shown  on  the  figure.  We  have  also 
computed  the  energy  for  4-particle  systems  for  a  variety  of 
lattice  sizes  and  find  agreement  with  Ref.  [25],  The  error  bands 
in  the  figure  provide  least-squares  estimates  for  the  one  a  error 
based  upon  quadratic  fits  to  the  finite-size  effects.  The  fits  are 
of  the  form  E /E-po  =  f  +  Ap  E3  +  Bp2/3.  For  the  interactions 
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FIG.  1 .  (Color  online)  The  calculated  ground-state  energy  shown 
as  the  value  of  f  vs  the  lattice  size  for  various  particle  numbers  and 
Hamiltonians. 


FIG.  2.  (Color  online)  The  ground-state  energy  as  a  function  of 
kFre :  comparison  of  DMC  and  AFQMC  results.  Dashed  lines  are 
DMC  results,  shifted  down  by  0.02  to  enable  comparison  of  the 
slopes. 


tuned  to  re  =  0,  a  fit  with  A  fixed  to  zero  is  used.  Including  a 
linear  coefficient  in  the  fit  yields  a  value  statistically  consistent 
with  zero. 

The  extrapolation  in  lattice  size  for  the  k2  and  Hubbard 
dispersions  show  opposite  slope  as  expected  from  the  opposite 
signs  of  their  effective  ranges  (see  Table  II).  The  extrapolation 
to  p  — >  0  is  consistent  with  £  =  0.372(0.005)  in  all  cases. 
Our  final  error  contains  a  statistical  component  and  the  errors 
associated  with  finite  population  sizes  and  finite  time-step 
errors.  This  value  is  below  previous  experiments  but  is  more 
compatible  with  recent  experimental  results  of  the  Zwierlein 
group  [8]. 

We  have  also  examined  the  behavior  of  the  energy  as 
a  function  of  kFre  for  finite  effective  ranges.  The  available 
length  scales  at  low  energy  restrict  £  to  be  a  function  of  kFa, 
kFre,  kpL ,  and  kpa.  At  unitarity  ( kFa)~l  —  0.  ( kFL)~ 3  is 
proportional  to  TV-1  so  for  a  sufficient  number  of  particles 
it  can  be  neglected.  The  lattice  has  inversion  symmetry  so 
we  expect  a  lowest  order  contribution  proportional  to  (kpa)2. 
Therefore,  the  only  linear  in  kp  correction  to  f  is  proportional 
to  kFre.  It  has  been  conjectured  [28,29]  that  the  slope  of  §  is 

TABLE  II.  Energy  extrapolations  to  infinite  volume,  zero  range 
limit  for  various  particle  numbers  N,  and  different  Hamiltonians.  The 
term  linear  in  the  effective  range,  A,  is  also  shown  where  it  is  not 
tuned  to  zero. 


Hamiltonian 

N 

? 

err 

A 

err 

(2) 

14 

0.39 

0.01 

0.21 

0.12 

38 

0.370 

0.005 

0.14 

0.04 

66 

0.374 

0.005 

0.11 

0.04 

(4) 

38 

0.372 

0.002 

48 

0.372 

0.003 

66 

0.372 

0.003 

Ah) 

ek 

4 

0.280 

0.004 

-0.28 

0.05 

38 

0.380 

0.005 

-0.17 

0.03 

48 

0.367 

0.005 

-0.05 

0.03 

66 

0.375 

0.005 

-0.13 

0.03 

universal  in  continuum  Hamiltonians:  £(re)  =  f  +  SkFre.  Of 
course,  a  finite-range  purely  attractive  interaction  is  subject 
to  collapse  for  a  many-particle  system,  but  a  small  repulsive 
many-body  interaction  or  the  lattice,  where  double  occupancy 
of  a  single  species  is  not  allowed,  is  enough  to  stabilize 
the  system.  Our  results  are  consistent  with  the  universality 
conjecture.  In  particular,  our  results  for  zero-effective-range 
approach  the  continuum  limit  with  a  slope  consistent  with 
zero. 

Figure  2  compares  the  AFQMC  results  for  the  ek  inter¬ 
action  with  the  DMC  results  [11,12]  for  various  values  of 
the  effective  range.  The  AFQMC  produces  somewhat  lower 
energies  than  the  DMC,  consistent  with  the  upper-bound  nature 
of  the  DMC  calculations.  For  the  slope  S  of  f  with  respect 
to  finite  re,  the  fit  to  the  N  =  66  AFQMC  results  yields 
S  =  0.11(0.03).  Similar  fits  to  the  AFQMC  data  with  the 
Hubbard  dispersion  ek  for  N  =  66  yield  a  linear  term  of 
S  =  0.12(0.03).  Both  are  in  agreement  with  the  DMC  results 
of  S  =  0.12(0.01)  [30]. 

In  summary,  we  have  developed  an  exact  lattice  technique 
to  treat  strongly  paired  fermion  systems.  We  find  that  the 
energy  of  the  Fermi  gas  at  unitarity  and  zero  effective  range 
is  §  =  0.372(0.005)  using  a  variety  of  interactions  tuned  to 
unitarity.  We  also  find  results  consistent  with  a  universal 
dependence  of  the  ground-state  energy  upon  the  effective 
range:  E/EFq  =  £  +  SkFre  +  •  •  •  with5  =  0.1 2(0.03)  for  the 
different  lattice  and  continuum  Hamiltonians.  The  method  we 
describe  should  be  useful  without  modification  for  the  entire 
BCS-BEC  transition  and  for  studying  many  other  properties 
of  cold  Fermi  gases.  It  can  also  be  applied  to  a  wide  variety  of 
problems  in  other  strongly  correlated  fermions  in  areas  ranging 
from  cold  atoms  to  condensed  matter  to  quantum  chemistry  to 
nuclear  physics. 
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